/*
 * Copyright 2009-2019 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
#define BOOST_TEST_MAIN

#define BOOST_TEST_MODULE orbitals_test
#include "votca/xtp/aomatrix.h"
#include <boost/test/unit_test.hpp>
#include <fstream>
#include <votca/xtp/convergenceacc.h>
#include <votca/xtp/orbitals.h>

using namespace votca::xtp;
using votca::Index;
BOOST_AUTO_TEST_SUITE(orbitals_test)
BOOST_AUTO_TEST_CASE(readxyztest) {
  std::ofstream xyzfile("molecule.xyz");
  xyzfile << " C 0.0 3 1" << std::endl;
  xyzfile << " methane" << std::endl;
  xyzfile << " C            .000000     .000000     .000000" << std::endl;
  xyzfile << " H            .629118     .629118     .629118" << std::endl;
  xyzfile << " H           -.629118    -.629118     .629118" << std::endl;
  xyzfile << " H            .629118    -.629118    -.629118" << std::endl;
  xyzfile << " H           -.629118     .629118    -.629118" << std::endl;
  xyzfile.close();

  Orbitals orb;

  BOOST_CHECK_THROW(orb.QMAtoms().LoadFromFile("molecule.xyz"),
                    std::runtime_error);
}

BOOST_AUTO_TEST_CASE(sortEnergies) {

  Orbitals orb;
  Eigen::VectorXd Energies = Eigen::VectorXd::LinSpaced(10, -5, 5);
  Eigen::VectorXd switched = Energies;
  switched(3) = Energies(5);
  switched(5) = Energies(3);
  orb.MOs().eigenvalues() = switched;
  orb.MOs().eigenvectors() = Eigen::MatrixXd::Zero(10, 10);
  orb.OrderMOsbyEnergy();
  bool issorted = Energies.isApprox(orb.MOs().eigenvalues(), 0.001);
  if (!issorted) {
    std::cout << "before" << std::endl;
    std::cout << Energies << std::endl;
    std::cout << "after" << std::endl;
    std::cout << orb.MOs().eigenvalues() << std::endl;
  }
  BOOST_CHECK_EQUAL(issorted, true);
}

BOOST_AUTO_TEST_CASE(densmat_test) {

  Orbitals orb;
  orb.setBasisSetSize(17);
  orb.setNumberOfOccupiedLevels(4);
  orb.setBSEindices(0, 9);
  orb.setNumberOfAlphaElectrons(5);
  orb.MOs().eigenvalues() = Eigen::VectorXd::Ones(17);
  orb.MOs().eigenvectors() = Eigen::MatrixXd::Zero(17, 17);
  orb.MOs().eigenvectors() << 0.9907442950821396, 0.20551857477099006,
      1.7023418853572578E-4, 9.074399061719398E-5, 3.837098252664792E-5,
      0.16974897759501614, -4.204899244661065E-5, 7.615775441246816E-5,
      4.676066062450439E-5, 2.884698428033532E-5, -4.481545968117419E-5,
      -1.3204180988296866E-4, -0.11680388320125862, 1.6189605485598983E-5,
      -5.371063226647223E-6, -1.2047015045305022E-4, 0.023240774591255682,
      0.08911047513015412, -0.20355821452675385, -2.4242411849295854E-4,
      -1.389669272938938E-4, -6.0728340934352946E-5, -0.07548903518938142,
      3.262334625557027E-5, -7.10564649816358E-5, -9.986644282246894E-5,
      1.6702170747853146E-4, -3.4464448598096653E-4, -0.0011237781000302065,
      0.16022562161339227, 2.5597538884079596E-4, -5.987456557611593E-4,
      0.0010289053314875896, 1.9979502865114958, 8.107959419149133E-6,
      -1.3746554617574666E-4, 0.30409982472519304, -0.04943981412502669,
      -0.18275903305786392, -1.4404163587992467E-4, 0.14257881572584336,
      -0.032034049310960475, 0.23956240448089597, 0.5387878318279876,
      0.07398471014538607, 0.4387442141244202, -0.0019929701691466604,
      -0.4478311747397298, 0.10346528783585213, 0.7800429786477285,
      4.413287003871695E-4, 1.1754110374944727E-5, -1.6786131622203274E-4,
      0.16187208368047987, 0.2472771093038189, 0.202328179665061,
      -3.820994854847884E-4, -0.15977453182924237, 0.19628611211192631,
      0.1213197329408795, -0.4388814545469389, 0.20281124109383114,
      0.5046453462735411, -0.002238807774688505, 0.5150776559618175,
      -0.6393382963354479, 0.38099226873671443, 6.6047695407522E-4,
      -5.048302992539679E-7, 1.3568022971562862E-5, -0.0982741302898521,
      0.2543655402180741, -0.23220810764458472, -1.944169891214517E-5,
      0.18144825026251918, 0.1980251300437136, -0.08147866661574833,
      0.07405511360833358, 0.6646693464629595, -0.2028284384893438,
      2.279342030052545E-4, -0.5940706293255381, -0.6323141290039537,
      -0.25763397034652824, -2.315661453010646E-5, -0.08143505806634609,
      -0.7045679157387013, -4.4787955662640034E-4, -2.0910866005311433E-4,
      -8.259360481466638E-5, -2.424093867046151, 6.780952279118335E-4,
      -0.0011263648445154907, -2.1613578613839522E-4, -3.6586704890798056E-4,
      3.4830029266647066E-4, 7.405837519085355E-4, 0.14778031305572506,
      -5.874372031528609E-4, 0.0012517213522246784, -0.0016624970087244567,
      -3.9307290475800984, -4.845180411930193E-5, -3.15232661550629E-4,
      0.29033356022603657, -0.047214119469099654, -0.1746229236165066,
      -0.0010623554092176837, 0.7157746323727515, -0.1610150829224024,
      1.207129610424005, -0.8861441619621693, -0.12082122795238953,
      -0.7253976412392744, 0.004702513046801041, 0.5445150195359842,
      -0.12577564486773046, -0.9479104485435895, -0.0029133750584969995,
      -6.97620867722761E-5, -4.518653975073802E-4, 0.1545537407856724,
      0.23619223591899627, 0.19332335260823175, -0.0022069712738471844,
      -0.8019796219162296, 0.9866731245664445, 0.6112535135056467,
      0.7197858303062914, -0.33414367829784547, -0.8313511890447793,
      0.005950677113076908, -0.6262908623920043, 0.7772100698368285,
      -0.46298826705698615, -0.004212137095734844, 3.025513825093562E-6,
      2.01122633970492E-5, -0.09381992712564609, 0.24295366644273125,
      -0.22187182348758724, -4.4849221276044244E-5, 0.9107435069506925,
      0.9953611683731745, -0.410462654769864, -0.11998173280053394,
      -1.0921777637934424, 0.3342503817907009, -4.2739718665356585E-4,
      0.7223406094651171, 0.7686726934386922, 0.313084531007334,
      1.7805853565722838E-4, 7.37158442829733E-4, -0.11947916702565467,
      0.0709877748234037, 0.19486216169399248, -0.19963587525346332,
      0.021970191718328234, -0.08852720779007611, -0.08752682338800259,
      -0.029183087590095056, 0.2404775527159037, 0.6601558927534981,
      0.12659420741971703, -0.667743162629014, 0.6085237963794983,
      0.5985008552792193, -0.20986989221784247, 0.15058182893557184,
      0.01608075222828349, 0.010011975370567302, 0.05470769366508042,
      0.15064492298573015, -0.15457287617332643, 0.9255720960213428,
      -1.206121064891608, -1.1996011487723761, -0.4056978894599689,
      -0.05298499024971882, -0.14303294036061362, -0.02485325238033535,
      0.2987439256438313, -1.013863139043371, -0.9971897529370147,
      0.3487761135053009, 0.7790198218470478, 7.138422212641049E-4,
      -0.1196094899026972, -0.2842166232364692, -0.04556245465364824,
      -0.007053279485955037, 0.02230987096947848, -0.0036852600571774035,
      0.023565132654508568, 0.12524632931555021, -0.17382906342060087,
      -0.060853394762126325, -0.6918002740597844, -0.6669226960233403,
      0.035864677061268974, -0.16556444794734262, 0.861747976464038,
      0.14795843313851584, 0.016063982425839392, 0.009859499054321354,
      -0.21897014622303052, -0.035202371814191405, -0.005465437420629127,
      0.9257581110182598, -0.04980711970695707, 0.32109962848034435,
      1.7221001218391445, 0.035905784702012784, 0.013498840982845527,
      0.14562342227832253, 0.3047901425448661, -0.05951259424946384,
      0.27553962984145547, -1.4351881046889507, 0.7788885658103283,
      7.409253682701713E-4, -0.11944350976971661, 0.07513719426852992,
      0.09761819951158109, 0.26015661269895035, 0.021951308529804293,
      0.11609692884630969, -0.043567065850759326, -0.03144638449326769,
      -0.5763047315656864, -0.12654072032163471, 0.40137211030812026,
      -0.6681794172539768, -0.7963825251648803, 0.30840285525063804,
      -0.2072016617793105, 0.15099996101012425, 0.016080755519097423,
      0.01002674093860076, 0.05792745695939718, 0.07547359788405249,
      0.20146343358941035, 0.9257451411073216, 1.5832800173788046,
      -0.5979818102877212, -0.4371140280860898, 0.12657732301631605,
      0.028622969925985273, -0.08503838333777329, 0.2984829177687489,
      1.3272620325313798, -0.5139436624743582, 0.3442485273058088,
      0.7788837381362927, 7.28527937995282E-4, -0.11955381155867925,
      0.13824691197079378, -0.24681524845383587, -0.05342022280053824,
      0.02203103650716836, -0.024102851640929047, 0.10805548656383365,
      -0.06369358394744931, 0.5104340186250341, -0.47455391846866973,
      0.15780550183151598, -0.6668890690900623, 0.15206929191456245,
      -0.740817337705803, -0.4483567934721803, 0.14962178851412328,
      0.016079430040486366, 0.00997356912486891, 0.10652865972788295,
      -0.19075159916388634, -0.04134322428836118, 0.9252734630049098,
      -0.3278217219514359, 1.4769403455748562, -0.8808246951973489,
      -0.1097870975651871, 0.10218284493845689, -0.03092608243612813,
      0.2996873524742791, -0.25325921205821206, 1.2339194146754904,
      0.7458659447098919, 0.7792524302192918;
  std::ofstream xyzfile("molecule.xyz");
  xyzfile << "5" << std::endl;
  xyzfile << " methane" << std::endl;
  xyzfile << " C            .000000     .000000     .000000" << std::endl;
  xyzfile << " H            .629118     .629118     .629118" << std::endl;
  xyzfile << " H           -.629118    -.629118     .629118" << std::endl;
  xyzfile << " H            .629118    -.629118    -.629118" << std::endl;
  xyzfile << " H           -.629118     .629118    -.629118" << std::endl;
  xyzfile.close();

  orb.QMAtoms().LoadFromFile("molecule.xyz");
  QMState s = QMState("n");
  Eigen::MatrixXd dmat_gs = orb.DensityMatrixFull(s);

  Eigen::MatrixXd dmat_gs_ref = Eigen::MatrixXd::Zero(17, 17);
  dmat_gs_ref << 2.04762, 0.0929013, 5.4126e-05, 5.42835e-05, 1.72817e-05,
      -0.450966, -0.000135298, -0.00022848, 2.64124e-05, -0.0475902, 0.0360251,
      -0.0478545, 0.0358023, -0.0475843, 0.0360186, -0.0476952, 0.0359624,
      0.0929013, 0.0987534, -7.62918e-05, -7.67761e-05, -2.86625e-05, 0.272328,
      -7.94404e-06, 3.09478e-05, -2.96854e-05, 0.0486847, -0.00127851,
      0.0489727, -0.00103507, 0.0486959, -0.00126519, 0.0488037, -0.00119335,
      5.4126e-05, -7.62918e-05, 0.189842, 0.0739999, -0.0849219, -5.93365e-05,
      0.181249, 0.0706451, -0.0810844, 0.0239397, 0.018375, -0.168322,
      -0.129699, 0.0360788, 0.0277662, 0.10852, 0.0836497, 5.42835e-05,
      -7.67761e-05, 0.0739999, 0.174697, 0.0939819, -1.3789e-05, 0.070644,
      0.166846, 0.0897801, 0.119392, 0.0922104, -0.114506, -0.0883027,
      0.0726428, 0.0560765, -0.0772667, -0.0598519, 1.72817e-05, -2.86625e-05,
      -0.0849219, 0.0939819, 0.148719, -3.73871e-05, -0.0810839, 0.0897811,
      0.142038, 0.0851767, 0.0658853, 0.03268, 0.0251299, 0.0348901, 0.0270105,
      -0.152738, -0.117979, -0.450966, 0.272328, -5.93365e-05, -1.3789e-05,
      -3.73871e-05, 1.0061, 0.000211774, 0.000410879, -4.64011e-05, 0.168097,
      -0.0168393, 0.168703, -0.0162988, 0.168083, -0.0168316, 0.168328,
      -0.0166886, -0.000135298, -7.94404e-06, 0.181249, 0.070644, -0.0810839,
      0.000211774, 0.173046, 0.0674414, -0.0774198, 0.022895, 0.017534,
      -0.160658, -0.123832, 0.034487, 0.0265019, 0.103657, 0.0798622,
      -0.00022848, 3.09478e-05, 0.0706451, 0.166846, 0.0897811, 0.000410879,
      0.0674414, 0.159348, 0.0857671, 0.114101, 0.0880616, -0.109268,
      -0.0843255, 0.0694466, 0.0535471, -0.0737506, -0.0571905, 2.64124e-05,
      -2.96854e-05, -0.0810844, 0.0897801, 0.142038, -4.64011e-05, -0.0774198,
      0.0857671, 0.135657, 0.08136, 0.0629346, 0.0311864, 0.0239829, 0.0333299,
      0.0258042, -0.145875, -0.112676, -0.0475902, 0.0486847, 0.0239397,
      0.119392, 0.0851767, 0.168097, 0.022895, 0.114101, 0.08136, 0.114573,
      0.0641084, -0.0295259, -0.04714, 0.0772549, 0.0352659, -0.0479928,
      -0.0615756, 0.0360251, -0.00127851, 0.018375, 0.0922104, 0.0658853,
      -0.0168393, 0.017534, 0.0880616, 0.0629346, 0.0641084, 0.0520913,
      -0.0471973, -0.0338508, 0.0352646, 0.0297955, -0.0616071, -0.0450988,
      -0.0478545, 0.0489727, -0.168322, -0.114506, 0.03268, 0.168703, -0.160658,
      -0.109268, 0.0311864, -0.0295259, -0.0471973, 0.194324, 0.125342,
      -0.0230317, -0.042181, -0.0274925, -0.0455351, 0.0358023, -0.00103507,
      -0.129699, -0.0883027, 0.0251299, -0.0162988, -0.123832, -0.0843255,
      0.0239829, -0.04714, -0.0338508, 0.125342, 0.0990848, -0.0421099,
      -0.0299681, -0.045501, -0.0325101, -0.0475843, 0.0486959, 0.0360788,
      0.0726428, 0.0348901, 0.168083, 0.034487, 0.0694466, 0.0333299, 0.0772549,
      0.0352646, -0.0230317, -0.0421099, 0.0588844, 0.0210688, 0.00114858,
      -0.0235919, 0.0360186, -0.00126519, 0.0277662, 0.0560765, 0.0270105,
      -0.0168316, 0.0265019, 0.0535471, 0.0258042, 0.0352659, 0.0297955,
      -0.042181, -0.0299681, 0.0210688, 0.018822, -0.0236135, -0.0157344,
      -0.0476952, 0.0488037, 0.10852, -0.0772667, -0.152738, 0.168328, 0.103657,
      -0.0737506, -0.145875, -0.0479928, -0.0616071, -0.0274925, -0.045501,
      0.00114858, -0.0236135, 0.188647, 0.121254, 0.0359624, -0.00119335,
      0.0836497, -0.0598519, -0.117979, -0.0166886, 0.0798622, -0.0571905,
      -0.112676, -0.0615756, -0.0450988, -0.0455351, -0.0325101, -0.0235919,
      -0.0157344, 0.121254, 0.0961851;

  bool check_dmat_gs = dmat_gs.isApprox(dmat_gs_ref, 0.0001);
  if (!check_dmat_gs) {
    std::cout << "Result gs" << std::endl;
    std::cout << dmat_gs << std::endl;
    std::cout << "Ref" << std::endl;
    std::cout << dmat_gs_ref << std::endl;
  }
  BOOST_CHECK_EQUAL(check_dmat_gs, 1);
  orb.setTDAApprox(false);
  orb.BSESinglets().eigenvectors() = Eigen::MatrixXd::Zero(25, 2);

  orb.BSESinglets().eigenvectors() << -9.172325499316926E-6,
      6.465589452839982E-6, 3.7017109343430703E-4, -9.603903429750687E-5,
      -1.1309760027198494E-4, -3.573759093655276E-4, -1.3654831774463808E-5,
      5.610635944952427E-5, 8.078036128870584E-4, 2.3508225792613987E-4,
      -4.1389689455982744E-4, 3.470240530817982E-4, 0.010935879545252295,
      -0.0029505107898146012, -0.003337063524453458, -0.010843128041228025,
      -6.022396289068329E-4, 0.0017252963787394244, 0.012406560343731091,
      0.0036243803634561714, -0.0534328659114822, 0.11790152418882739,
      -0.014341255912595252, 0.002244216527436879, 0.006333920883336163,
      0.022446822466352594, -0.013418752450867943, -0.005038670208008163,
      -0.012170499758537307, -0.006597144566176471, -0.2818369874026935,
      -0.9534596388889569, -0.020055422721015168, 0.019417142136606193,
      -0.019391880882433613, -0.009364499588916334, 0.008461551487656474,
      0.022977612625792006, -0.0036159338191032687, 0.01947317587720039,
      -0.9569542694576304, 0.2745425547752314, -0.025652859518621507,
      -0.019427980787681504, 0.02011572854778333, -0.022825383620618405,
      0.016907164519736315, -8.594057567669325E-4, -0.01566776126491342,
      -0.0036991931035359783;

  orb.BSESinglets().eigenvectors2() = Eigen::MatrixXd::Zero(25, 2);
  orb.BSESinglets().eigenvectors2() << 8.507035131871254E-6,
      -7.771078153320694E-6, -6.592414632166431E-4, 1.8308223920734227E-4,
      1.9536776547437147E-4, 6.525900973931339E-4, 3.278603510053265E-5,
      -9.958063857447739E-5, -0.0011044648385455508, -3.198831712590277E-4,
      1.3986989346956674E-4, -1.0389121949776523E-4, 0.011288355421213793,
      -0.003187756508136955, -0.0033091549363569992, -0.011220902898107834,
      -5.812495659517037E-4, 0.0016919649089329493, -7.785385054188999E-4,
      -2.135292962729838E-4, -0.0014926684696797512, 0.0032197968277020163,
      -0.007005339684942067, 0.001062499099375741, 0.0030075906847693325,
      0.00938345421251267, -0.0032541108414566823, -0.0036930063078656388,
      0.0017512723306478816, 9.058993399572328E-4, -0.007678356616119832,
      -0.02591115391341041, -0.009034931850146257, 0.01014588913521307,
      -0.007764559719931153, -0.004407761436169589, 0.0031656406702596303,
      0.009697260357617128, 4.730611502982343E-4, -0.0026997355242154015,
      -0.026013063963669262, 0.007447618201187667, -0.011436219477735475,
      -0.00784454627778686, 0.009258965011894391, -0.009898606114539346,
      0.006465113614564366, -9.868851256296352E-4, 0.0021031504117920274,
      5.800186595669497E-4;

  QMState s1 = QMState("s1");
  Eigen::MatrixXd dmat_s1 = orb.DensityMatrixFull(s1);

  Eigen::MatrixXd dmat_s1_ref = Eigen::MatrixXd::Zero(17, 17);
  dmat_s1_ref << 2.0474936814, 0.093039309013, 0.00027930870559,
      0.00013412117551, 0.0002236273963, -0.45057423895, 0.00032447883031,
      9.6414976839e-05, 0.0011453614901, -0.047569979558, 0.034564826654,
      -0.047863976862, 0.03605787982, -0.047474712999, 0.036553682373,
      -0.047481413274, 0.036747708877, 0.093039309013, 0.098607081953,
      -0.00023279635112, -0.00016344664426, -0.00011627503089, 0.27186176359,
      -0.00025977030654, -0.00017164932468, -0.00051684983536, 0.048597963266,
      -0.00065651947538, 0.048987165838, -0.001100443728, 0.048577467033,
      -0.0014975538044, 0.048647287156, -0.0015531572202, 0.00027930870559,
      -0.00023279635112, 0.19655470027, 0.062483612277, -0.049653026698,
      -0.0020243015878, 0.25644278848, -0.014767477527, 0.044645570493,
      0.019014899921, -0.12447218624, -0.16362165863, -0.13113328777,
      0.052109906179, 0.21355547661, 0.092493654522, 0.044658112006,
      0.00013412117551, -0.00016344664426, 0.062483612277, 0.13805903626,
      0.01374586465, -0.00059018682089, -0.022483358425, 0.2342384044,
      -0.082817359094, 0.085120139552, 0.20997986792, -0.098492150743,
      -0.060177827507, 0.030698425252, -0.18697362472, -0.017306020655,
      0.037565194415, 0.0002236273963, -0.00011627503089, -0.049653026698,
      0.01374586465, 0.12583059032, -0.003453325581, 0.050907693873,
      -0.086367571801, 0.26723896845, 0.022964754867, -0.20191414976,
      0.042207782543, 0.027471211902, 0.032543541477, 0.27453629016,
      -0.097659837212, -0.094773081646, -0.45057423895, 0.27186176359,
      -0.0020243015878, -0.00059018682089, -0.003453325581, 1.0056879047,
      -0.0051247721656, -0.0037258816235, -0.016431185902, 0.16890751016,
      0.0040432062202, 0.16859040691, -0.020894795359, 0.16767752387,
      -0.024247448135, 0.16679382414, -0.027502187342, 0.00032447883031,
      -0.00025977030654, 0.25644278848, -0.022483358425, 0.050907693873,
      -0.0051247721656, 0.62174241904, -0.45110587693, 0.53643247096,
      -0.030214912981, -0.75698204974, -0.1575608267, -0.16710662043,
      0.11609541262, 1.0468688747, 0.07186019552, -0.11470344006,
      9.6414976839e-05, -0.00017164932468, -0.014767477527, 0.2342384044,
      -0.086367571801, -0.0037258816235, -0.45110587693, 0.7827128442,
      -0.58341974054, 0.13254479911, 0.85226818291, -0.090846092114,
      0.025514492598, -0.051804518261, -1.166193648, 0.010706772158,
      0.29377112043, 0.0011453614901, -0.00051684983536, 0.044645570493,
      -0.082817359094, 0.26723896845, -0.016431185902, 0.53643247096,
      -0.58341974054, 0.9908033055, -0.048178244442, -1.1351900819,
      0.039805445697, 0.0015876844895, 0.11327220026, 1.3415541164,
      -0.1044896076, -0.18250168934, -0.047569979558, 0.048597963266,
      0.019014899921, 0.085120139552, 0.022964754867, 0.16890751016,
      -0.030214912981, 0.13254479911, -0.048178244442, 0.087000185799,
      0.14743799655, -0.018593966534, -0.034558271129, 0.048726392738,
      -0.10992334557, -0.0032313074277, -0.014314664106, 0.034564826654,
      -0.00065651947538, -0.12447218624, 0.20997986792, -0.20191414976,
      0.0040432062202, -0.75698204974, 0.85226818291, -1.1351900819,
      0.14743799655, 1.5258086619, -0.038222860143, 0.0082651239446,
      -0.10916739961, -1.6985408041, -0.010179663354, 0.13479752897,
      -0.047863976862, 0.048987165838, -0.16362165863, -0.098492150743,
      0.042207782543, 0.16859040691, -0.1575608267, -0.090846092114,
      0.039805445697, -0.018593966534, -0.038222860143, 0.18556977382,
      0.11990183188, -0.017133317747, -0.041766372367, -0.035473121585,
      -0.049165724862, 0.03605787982, -0.001100443728, -0.13113328777,
      -0.060177827507, 0.027471211902, -0.020894795359, -0.16710662043,
      0.025514492598, 0.0015876844895, -0.034558271129, 0.0082651239446,
      0.11990183188, 0.12690462965, -0.050224098204, -0.15534424247,
      -0.044067095492, 0.030332887428, -0.047474712999, 0.048577467033,
      0.052109906179, 0.030698425252, 0.032543541477, 0.16767752387,
      0.11609541262, -0.051804518261, 0.11327220026, 0.048726392738,
      -0.10916739961, -0.017133317747, -0.050224098204, 0.064125611198,
      0.19590547187, 0.018181459037, -0.045848940579, 0.036553682373,
      -0.0014975538044, 0.21355547661, -0.18697362472, 0.27453629016,
      -0.024247448135, 1.0468688747, -1.166193648, 1.3415541164, -0.10992334557,
      -1.6985408041, -0.041766372367, -0.15534424247, 0.19590547187,
      2.3098850938, -0.053821599569, -0.44099062588, -0.047481413274,
      0.048647287156, 0.092493654522, -0.017306020655, -0.097659837212,
      0.16679382414, 0.07186019552, 0.010706772158, -0.1044896076,
      -0.0032313074277, -0.010179663354, -0.035473121585, -0.044067095492,
      0.018181459037, -0.053821599569, 0.13469507677, 0.10041313642,
      0.036747708877, -0.0015531572202, 0.044658112006, 0.037565194415,
      -0.094773081646, -0.027502187342, -0.11470344006, 0.29377112043,
      -0.18250168934, -0.014314664106, 0.13479752897, -0.049165724862,
      0.030332887428, -0.045848940579, -0.44099062588, 0.10041313642,
      0.29490774526;

  bool check_dmat_s1 = dmat_s1.isApprox(dmat_s1_ref, 0.0001);
  if (!check_dmat_s1) {
    std::cout << "Result s1" << std::endl;
    std::cout << dmat_s1 << std::endl;
    std::cout << "Ref" << std::endl;
    std::cout << dmat_s1_ref << std::endl;
  }
  BOOST_CHECK_EQUAL(check_dmat_s1, 1);

  QMState n2s1 = QMState("n2s1");
  Eigen::MatrixXd dmat_n2s1 = orb.DensityMatrixFull(n2s1);
  Eigen::MatrixXd dmat_n2s1_ref = Eigen::MatrixXd::Zero(17, 17);
  dmat_n2s1_ref << -0.000397846, 0.000176625, -0.00881468, 0.00292497,
      -5.88628e-05, 0.00567658, -0.0158532, -0.00993011, 0.00907284,
      -0.00270861, 0.000263923, -0.000924017, -0.0295645, 0.00575722,
      0.00290312, -0.00234934, 0.0177591, 0.000320067, -0.000142077, 0.00842665,
      -0.00278916, 0.000112243, -0.00456611, 0.0156573, 0.00916967, -0.00851408,
      0.00256481, -0.00076839, 0.000922021, 0.0281029, -0.00544032, -0.00286575,
      0.0021369, -0.0175249, -0.00119257, 0.000534394, 0.00182468, -0.0294704,
      -0.0130282, 0.0170345, 0.0592909, -0.18966, -0.0559797, 0.00141623,
      0.0684471, -0.00279673, -0.0743399, 0.0215782, 0.173514, -0.0209129,
      -0.19341, -0.000825211, 0.000364902, -0.046019, 0.0405825, -0.0800878,
      0.0117652, -0.23343, 0.204235, -0.397217, 0.0365988, 0.505462,
      -0.000632154, -0.0064756, -0.0364718, -0.50527, 0.000110369, -0.0117897,
      0.000219408, -0.000102567, -0.0415354, 0.0585534, -0.0595295, -0.00315313,
      -0.250057, 0.327208, -0.302048, 0.0307742, 0.386554, 0.0016499, 0.0528536,
      -0.0487586, -0.576802, 0.0165534, 0.141939, 0.00113314, -0.000503003,
      0.0292259, -0.00966418, 0.000269493, -0.0161657, 0.053978, 0.0321033,
      -0.0301292, 0.00896316, -0.00185594, 0.0031865, 0.09765, -0.0189548,
      -0.010399, 0.00745341, -0.060809, -0.00113828, 0.000510066, 0.00175184,
      -0.0281414, -0.0124351, 0.016259, 0.0566329, -0.18108, -0.0534369,
      0.00135285, 0.065327, -0.00266934, -0.0709506, 0.0205986, 0.165684,
      -0.0199648, -0.184674, -0.000787497, 0.000348224, -0.0439426, 0.0387632,
      -0.0764923, 0.0112275, -0.222946, 0.19512, -0.379399, 0.0349603, 0.482771,
      -0.00060196, -0.00613095, -0.0348473, -0.482638, 0.000112355, -0.0112504,
      0.000209371, -9.78784e-05, -0.0396733, 0.0559237, -0.0568627, -0.00300889,
      -0.238837, 0.312508, -0.288514, 0.0293946, 0.369235, 0.00157545,
      0.0504712, -0.046569, -0.550915, 0.0158079, 0.135546, -0.000231385,
      0.000100807, -0.0300844, 0.0339378, -0.0581988, 0.00329184, -0.176959,
      0.191403, -0.295432, 0.0288963, 0.369598, 0.000532762, 0.0241301,
      -0.034354, -0.412182, 0.0048565, 0.0132905, -0.000342812, 0.000150699,
      -0.0274952, 0.0276502, -0.045003, 0.00488358, -0.144628, 0.143375,
      -0.223878, 0.0210264, 0.285828, -4.91689e-05, 0.00451078, -0.0238081,
      -0.317125, 0.00268425, 0.0192386, 0.00136974, -0.000611345, 0.0190123,
      0.00703061, 0.0367438, -0.0195544, 0.0426127, 0.0803459, 0.170366,
      -0.0118181, -0.221516, 0.00286778, 0.0749488, -0.0073028, 0.03646,
      0.0170358, 0.139817, 0.000891328, -0.000398214, 0.0104389, 0.00679366,
      0.02831, -0.0127262, 0.025138, 0.0571348, 0.135819, -0.0104218, -0.170651,
      0.00174826, 0.0435996, -0.00286107, 0.0298649, 0.0120438, 0.116519,
      -0.000180458, 7.97071e-05, -0.0134737, 0.0136744, -0.032556, 0.00257215,
      -0.0826179, 0.0810933, -0.166509, 0.0162635, 0.204982, 0.00019225,
      0.0114336, -0.017132, -0.19887, 0.000602948, -0.0215253, -0.000303473,
      0.000134402, -0.0146546, 0.0119854, -0.0251818, 0.00432765, -0.0716987,
      0.0580983, -0.124224, 0.0112615, 0.158582, -0.000312438, -0.00530481,
      -0.010495, -0.152229, -0.000604447, -0.00768288, -0.000195069,
      9.22122e-05, 0.044308, -0.0611731, 0.0540908, 0.00280769, 0.253367,
      -0.331024, 0.270639, -0.0272218, -0.353694, -0.00143346, -0.0443903,
      0.0459175, 0.567096, -0.0174666, -0.173028, -0.000314098, 0.000143768,
      0.0300096, -0.0458553, 0.0418061, 0.00450015, 0.187987, -0.260365,
      0.213736, -0.022355, -0.273329, -0.00156761, -0.0484174, 0.0382315,
      0.439776, -0.0145603, -0.124678;

  bool check_dmat_n2s1 = dmat_n2s1.isApprox(dmat_n2s1_ref, 0.0001);
  if (!check_dmat_n2s1) {
    std::cout << "Result n2s1" << std::endl;
    std::cout << dmat_n2s1 << std::endl;
    std::cout << "Ref" << std::endl;
    std::cout << dmat_n2s1_ref << std::endl;
  }

  BOOST_CHECK_EQUAL(check_dmat_n2s1, 1);
}

BOOST_AUTO_TEST_CASE(dipole_test) {

  std::ofstream xyzfile("molecule.xyz");
  xyzfile << " 5" << std::endl;
  xyzfile << " methane" << std::endl;
  xyzfile << " C            .000000     .000000     .000000" << std::endl;
  xyzfile << " H            .629118     .629118     .629118" << std::endl;
  xyzfile << " H           -.629118    -.629118     .629118" << std::endl;
  xyzfile << " H            .629118    -.629118    -.629118" << std::endl;
  xyzfile << " H           -.629118     .629118    -.629118" << std::endl;
  xyzfile.close();

  std::ofstream basisfile("3-21G.xml");
  basisfile << "<basis name=\"3-21G\">" << std::endl;
  basisfile << "  <element name=\"H\">" << std::endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << std::endl;
  basisfile << "      <constant decay=\"5.447178e+00\">" << std::endl;
  basisfile << "        <contractions factor=\"1.562850e-01\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "      <constant decay=\"8.245470e-01\">" << std::endl;
  basisfile << "        <contractions factor=\"9.046910e-01\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "    </shell>" << std::endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << std::endl;
  basisfile << "      <constant decay=\"1.831920e-01\">" << std::endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "    </shell>" << std::endl;
  basisfile << "  </element>" << std::endl;
  basisfile << "  <element name=\"C\">" << std::endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << std::endl;
  basisfile << "      <constant decay=\"1.722560e+02\">" << std::endl;
  basisfile << "        <contractions factor=\"6.176690e-02\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "      <constant decay=\"2.591090e+01\">" << std::endl;
  basisfile << "        <contractions factor=\"3.587940e-01\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "      <constant decay=\"5.533350e+00\">" << std::endl;
  basisfile << "        <contractions factor=\"7.007130e-01\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "    </shell>" << std::endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << std::endl;
  basisfile << "      <constant decay=\"3.664980e+00\">" << std::endl;
  basisfile << "        <contractions factor=\"-3.958970e-01\" type=\"S\"/>"
            << std::endl;
  basisfile << "        <contractions factor=\"2.364600e-01\" type=\"P\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "      <constant decay=\"7.705450e-01\">" << std::endl;
  basisfile << "        <contractions factor=\"1.215840e+00\" type=\"S\"/>"
            << std::endl;
  basisfile << "        <contractions factor=\"8.606190e-01\" type=\"P\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "    </shell>" << std::endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << std::endl;
  basisfile << "      <constant decay=\"1.958570e-01\">" << std::endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << std::endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"P\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "    </shell>" << std::endl;
  basisfile << "  </element>" << std::endl;
  basisfile << "</basis>" << std::endl;
  basisfile.close();

  Orbitals orbitals;
  orbitals.QMAtoms().LoadFromFile("molecule.xyz");
  BasisSet basis;
  basis.Load("3-21G.xml");
  orbitals.setDFTbasisName("3-21G.xml");
  AOBasis aobasis;
  aobasis.Fill(basis, orbitals.QMAtoms());

  orbitals.setBasisSetSize(17);
  orbitals.setNumberOfOccupiedLevels(4);
  Eigen::MatrixXd& MOs = orbitals.MOs().eigenvectors();
  orbitals.MOs().eigenvalues() = Eigen::VectorXd::Ones(17);
  MOs = Eigen::MatrixXd::Zero(17, 17);
  MOs << -0.00761992, -4.69664e-13, 8.35009e-15, -1.15214e-14, -0.0156169,
      -2.23157e-12, 1.52916e-14, 2.10997e-15, 8.21478e-15, 3.18517e-15,
      2.89043e-13, -0.00949189, 1.95787e-12, 1.22168e-14, -2.63092e-15,
      -0.22227, 1.00844, 0.233602, -3.18103e-12, 4.05093e-14, -4.70943e-14,
      0.1578, 4.75897e-11, -1.87447e-13, -1.02418e-14, 6.44484e-14, -2.6602e-14,
      6.5906e-12, -0.281033, -6.67755e-12, 2.70339e-14, -9.78783e-14, -1.94373,
      -0.36629, -1.63678e-13, -0.22745, -0.054851, 0.30351, 3.78688e-11,
      -0.201627, -0.158318, -0.233561, -0.0509347, -0.650424, 0.452606,
      -5.88565e-11, 0.453936, -0.165715, -0.619056, 7.0149e-12, 2.395e-14,
      -4.51653e-14, -0.216509, 0.296975, -0.108582, 3.79159e-11, -0.199301,
      0.283114, -0.0198557, 0.584622, 0.275311, 0.461431, -5.93732e-11,
      0.453057, 0.619523, 0.166374, 7.13235e-12, 2.56811e-14, -9.0903e-14,
      -0.21966, -0.235919, -0.207249, 3.75979e-11, -0.199736, -0.122681,
      0.255585, -0.534902, 0.362837, 0.461224, -5.91028e-11, 0.453245,
      -0.453298, 0.453695, 7.01644e-12, 2.60987e-14, 0.480866, 1.8992e-11,
      -2.56795e-13, 4.14571e-13, 2.2709, 4.78615e-10, -2.39153e-12,
      -2.53852e-13, -2.15605e-13, -2.80359e-13, 7.00137e-12, 0.145171,
      -1.96136e-11, -2.24876e-13, -2.57294e-14, 4.04176, 0.193617, -1.64421e-12,
      -0.182159, -0.0439288, 0.243073, 1.80753e-10, -0.764779, -0.600505,
      -0.885907, 0.0862014, 1.10077, -0.765985, 6.65828e-11, -0.579266,
      0.211468, 0.789976, -1.41532e-11, -1.29659e-13, -1.64105e-12, -0.173397,
      0.23784, -0.0869607, 1.80537e-10, -0.755957, 1.07386, -0.0753135,
      -0.989408, -0.465933, -0.78092, 6.72256e-11, -0.578145, -0.790571,
      -0.212309, -1.42443e-11, -1.31306e-13, -1.63849e-12, -0.17592, -0.188941,
      -0.165981, 1.79403e-10, -0.757606, -0.465334, 0.969444, 0.905262,
      -0.61406, -0.78057, 6.69453e-11, -0.578385, 0.578453, -0.578959,
      -1.40917e-11, -1.31002e-13, 0.129798, -0.274485, 0.00256652, -0.00509635,
      -0.0118465, 0.141392, -0.000497905, -0.000510338, -0.000526798,
      -0.00532572, 0.596595, 0.65313, -0.964582, -0.000361559, -0.000717866,
      -0.195084, 0.0246232, 0.0541331, -0.255228, 0.00238646, -0.0047388,
      -0.88576, 1.68364, -0.00592888, -0.00607692, -9.5047e-05, -0.000960887,
      0.10764, -0.362701, 1.53456, 0.000575205, 0.00114206, -0.793844,
      -0.035336, 0.129798, 0.0863299, -0.0479412, 0.25617, -0.0118465,
      -0.0464689, 0.0750316, 0.110468, -0.0436647, -0.558989, -0.203909,
      0.65313, 0.320785, 0.235387, 0.878697, -0.195084, 0.0246232, 0.0541331,
      0.0802732, -0.0445777, 0.238198, -0.88576, -0.553335, 0.893449, 1.31541,
      -0.00787816, -0.100855, -0.0367902, -0.362701, -0.510338, -0.374479,
      -1.39792, -0.793844, -0.035336, 0.129798, 0.0927742, -0.197727, -0.166347,
      -0.0118465, -0.0473592, 0.0582544, -0.119815, -0.463559, 0.320126,
      -0.196433, 0.65313, 0.321765, 0.643254, -0.642737, -0.195084, 0.0246232,
      0.0541331, 0.0862654, -0.183855, -0.154677, -0.88576, -0.563936, 0.693672,
      -1.42672, -0.0836372, 0.0577585, -0.0354411, -0.362701, -0.511897,
      -1.02335, 1.02253, -0.793844, -0.035336, 0.129798, 0.0953806, 0.243102,
      -0.0847266, -0.0118465, -0.0475639, -0.132788, 0.00985812, 0.507751,
      0.244188, -0.196253, 0.65313, 0.322032, -0.87828, -0.235242, -0.195084,
      0.0246232, 0.0541331, 0.088689, 0.226046, -0.0787824, -0.88576, -0.566373,
      -1.58119, 0.117387, 0.0916104, 0.0440574, -0.0354087, -0.362701,
      -0.512321, 1.39726, 0.374248, -0.793844, -0.035336;

  orbitals.setBSEindices(0, 16);
  orbitals.setTDAApprox(true);

  Eigen::MatrixXd spsi_ref = Eigen::MatrixXd::Zero(60, 1);
  spsi_ref << -0.000150849, 0.00516987, 0.0511522, 0.00428958, -0.00966668,
      -0.000155227, 1.02978e-08, 5.82225e-05, -0.00216177, 0.00907102,
      6.297e-09, -9.84993e-11, 0.00159727, 0.0039042, 0.0481196, 0.00495382,
      -0.0106013, 0.00025141, -0.000155626, -0.000382828, -0.00322057,
      0.0124251, 1.32177e-05, 6.794e-07, -0.0153713, 0.0200649, -0.067081,
      -0.0122678, 0.0117612, -0.00358901, 0.00605007, 0.00404793, 0.0108884,
      -0.0151075, -0.000513827, -2.64139e-05, -0.0466653, 0.0672016, 0.021747,
      -0.0115096, -0.0124868, -0.0115055, 0.0187191, 0.0124754, 0.0149534,
      0.0112807, -0.00158977, -8.17254e-05, -0.00290157, 0.0994541, 0.984029,
      0.017835, -0.0401912, -0.000645537, -7.54896e-08, -5.91055e-05,
      0.00219348, -0.00920484, 1.82832e-08, 5.56223e-11;

  orbitals.BSESinglets().eigenvectors() = spsi_ref;
  QMState state_trans = QMState("n2s1");

  Eigen::Vector3d res_trans = orbitals.CalcElDipole(state_trans);
  Eigen::Vector3d ref_trans = Eigen::Vector3d::Zero();
  ref_trans << 0.118565, 0.0444239, -0.0505149;

  bool check_trans = ref_trans.isApprox(res_trans, 0.0001);
  if (!check_trans) {
    std::cout << "Result transition dipole" << std::endl;
    std::cout << res_trans << std::endl;
    std::cout << "Ref transition dipole" << std::endl;
    std::cout << ref_trans << std::endl;
  }
  BOOST_CHECK_EQUAL(check_trans, 1);

  QMState state_s1 = QMState("s1");
  Eigen::Vector3d res_s1 = orbitals.CalcElDipole(state_s1);
  Eigen::Vector3d ref_s1 = Eigen::Vector3d::Zero();
  ref_s1 << -0.15153501734, -0.42406579479, 0.033954362839;
  bool check_s1 = ref_s1.isApprox(res_s1, 0.0001);
  if (!check_s1) {
    std::cout << "Result s1 dipole" << std::endl;
    std::cout << res_s1 << std::endl;
    std::cout << "Ref s1 dipole" << std::endl;
    std::cout << ref_s1 << std::endl;
  }
  BOOST_CHECK_EQUAL(check_s1, 1);
}

BOOST_AUTO_TEST_CASE(osc_strength) {

  std::ofstream xyzfile("molecule.xyz");
  xyzfile << " 5" << std::endl;
  xyzfile << " methane" << std::endl;
  xyzfile << " C            .000000     .000000     .000000" << std::endl;
  xyzfile << " H            .629118     .629118     .629118" << std::endl;
  xyzfile << " H           -.629118    -.629118     .629118" << std::endl;
  xyzfile << " H            .629118    -.629118    -.629118" << std::endl;
  xyzfile << " H           -.629118     .629118    -.629118" << std::endl;
  xyzfile.close();

  std::ofstream basisfile("3-21G.xml");
  basisfile << "<basis name=\"3-21G\">" << std::endl;
  basisfile << "  <element name=\"H\">" << std::endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << std::endl;
  basisfile << "      <constant decay=\"5.447178e+00\">" << std::endl;
  basisfile << "        <contractions factor=\"1.562850e-01\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "      <constant decay=\"8.245470e-01\">" << std::endl;
  basisfile << "        <contractions factor=\"9.046910e-01\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "    </shell>" << std::endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << std::endl;
  basisfile << "      <constant decay=\"1.831920e-01\">" << std::endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "    </shell>" << std::endl;
  basisfile << "  </element>" << std::endl;
  basisfile << "  <element name=\"C\">" << std::endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << std::endl;
  basisfile << "      <constant decay=\"1.722560e+02\">" << std::endl;
  basisfile << "        <contractions factor=\"6.176690e-02\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "      <constant decay=\"2.591090e+01\">" << std::endl;
  basisfile << "        <contractions factor=\"3.587940e-01\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "      <constant decay=\"5.533350e+00\">" << std::endl;
  basisfile << "        <contractions factor=\"7.007130e-01\" type=\"S\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "    </shell>" << std::endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << std::endl;
  basisfile << "      <constant decay=\"3.664980e+00\">" << std::endl;
  basisfile << "        <contractions factor=\"-3.958970e-01\" type=\"S\"/>"
            << std::endl;
  basisfile << "        <contractions factor=\"2.364600e-01\" type=\"P\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "      <constant decay=\"7.705450e-01\">" << std::endl;
  basisfile << "        <contractions factor=\"1.215840e+00\" type=\"S\"/>"
            << std::endl;
  basisfile << "        <contractions factor=\"8.606190e-01\" type=\"P\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "    </shell>" << std::endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << std::endl;
  basisfile << "      <constant decay=\"1.958570e-01\">" << std::endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << std::endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"P\"/>"
            << std::endl;
  basisfile << "      </constant>" << std::endl;
  basisfile << "    </shell>" << std::endl;
  basisfile << "  </element>" << std::endl;
  basisfile << "</basis>" << std::endl;
  basisfile.close();

  Orbitals orb;
  orb.QMAtoms().LoadFromFile("molecule.xyz");
  orb.setDFTbasisName("3-21G.xml");

  QMState s("s1");
  orb.setBasisSetSize(17);
  orb.setNumberOfOccupiedLevels(4);
  orb.setBSEindices(0, 16);
  orb.setTDAApprox(true);

  Eigen::MatrixXd& MOs = orb.MOs().eigenvectors();
  orb.MOs().eigenvalues() = Eigen::VectorXd::Ones(17);
  MOs = Eigen::MatrixXd::Zero(17, 17);
  MOs << -0.00761992, -4.69664e-13, 8.35009e-15, -1.15214e-14, -0.0156169,
      -2.23157e-12, 1.52916e-14, 2.10997e-15, 8.21478e-15, 3.18517e-15,
      2.89043e-13, -0.00949189, 1.95787e-12, 1.22168e-14, -2.63092e-15,
      -0.22227, 1.00844, 0.233602, -3.18103e-12, 4.05093e-14, -4.70943e-14,
      0.1578, 4.75897e-11, -1.87447e-13, -1.02418e-14, 6.44484e-14, -2.6602e-14,
      6.5906e-12, -0.281033, -6.67755e-12, 2.70339e-14, -9.78783e-14, -1.94373,
      -0.36629, -1.63678e-13, -0.22745, -0.054851, 0.30351, 3.78688e-11,
      -0.201627, -0.158318, -0.233561, -0.0509347, -0.650424, 0.452606,
      -5.88565e-11, 0.453936, -0.165715, -0.619056, 7.0149e-12, 2.395e-14,
      -4.51653e-14, -0.216509, 0.296975, -0.108582, 3.79159e-11, -0.199301,
      0.283114, -0.0198557, 0.584622, 0.275311, 0.461431, -5.93732e-11,
      0.453057, 0.619523, 0.166374, 7.13235e-12, 2.56811e-14, -9.0903e-14,
      -0.21966, -0.235919, -0.207249, 3.75979e-11, -0.199736, -0.122681,
      0.255585, -0.534902, 0.362837, 0.461224, -5.91028e-11, 0.453245,
      -0.453298, 0.453695, 7.01644e-12, 2.60987e-14, 0.480866, 1.8992e-11,
      -2.56795e-13, 4.14571e-13, 2.2709, 4.78615e-10, -2.39153e-12,
      -2.53852e-13, -2.15605e-13, -2.80359e-13, 7.00137e-12, 0.145171,
      -1.96136e-11, -2.24876e-13, -2.57294e-14, 4.04176, 0.193617, -1.64421e-12,
      -0.182159, -0.0439288, 0.243073, 1.80753e-10, -0.764779, -0.600505,
      -0.885907, 0.0862014, 1.10077, -0.765985, 6.65828e-11, -0.579266,
      0.211468, 0.789976, -1.41532e-11, -1.29659e-13, -1.64105e-12, -0.173397,
      0.23784, -0.0869607, 1.80537e-10, -0.755957, 1.07386, -0.0753135,
      -0.989408, -0.465933, -0.78092, 6.72256e-11, -0.578145, -0.790571,
      -0.212309, -1.42443e-11, -1.31306e-13, -1.63849e-12, -0.17592, -0.188941,
      -0.165981, 1.79403e-10, -0.757606, -0.465334, 0.969444, 0.905262,
      -0.61406, -0.78057, 6.69453e-11, -0.578385, 0.578453, -0.578959,
      -1.40917e-11, -1.31002e-13, 0.129798, -0.274485, 0.00256652, -0.00509635,
      -0.0118465, 0.141392, -0.000497905, -0.000510338, -0.000526798,
      -0.00532572, 0.596595, 0.65313, -0.964582, -0.000361559, -0.000717866,
      -0.195084, 0.0246232, 0.0541331, -0.255228, 0.00238646, -0.0047388,
      -0.88576, 1.68364, -0.00592888, -0.00607692, -9.5047e-05, -0.000960887,
      0.10764, -0.362701, 1.53456, 0.000575205, 0.00114206, -0.793844,
      -0.035336, 0.129798, 0.0863299, -0.0479412, 0.25617, -0.0118465,
      -0.0464689, 0.0750316, 0.110468, -0.0436647, -0.558989, -0.203909,
      0.65313, 0.320785, 0.235387, 0.878697, -0.195084, 0.0246232, 0.0541331,
      0.0802732, -0.0445777, 0.238198, -0.88576, -0.553335, 0.893449, 1.31541,
      -0.00787816, -0.100855, -0.0367902, -0.362701, -0.510338, -0.374479,
      -1.39792, -0.793844, -0.035336, 0.129798, 0.0927742, -0.197727, -0.166347,
      -0.0118465, -0.0473592, 0.0582544, -0.119815, -0.463559, 0.320126,
      -0.196433, 0.65313, 0.321765, 0.643254, -0.642737, -0.195084, 0.0246232,
      0.0541331, 0.0862654, -0.183855, -0.154677, -0.88576, -0.563936, 0.693672,
      -1.42672, -0.0836372, 0.0577585, -0.0354411, -0.362701, -0.511897,
      -1.02335, 1.02253, -0.793844, -0.035336, 0.129798, 0.0953806, 0.243102,
      -0.0847266, -0.0118465, -0.0475639, -0.132788, 0.00985812, 0.507751,
      0.244188, -0.196253, 0.65313, 0.322032, -0.87828, -0.235242, -0.195084,
      0.0246232, 0.0541331, 0.088689, 0.226046, -0.0787824, -0.88576, -0.566373,
      -1.58119, 0.117387, 0.0916104, 0.0440574, -0.0354087, -0.362701,
      -0.512321, 1.39726, 0.374248, -0.793844, -0.035336;

  Eigen::VectorXd se_ref = Eigen::VectorXd::Zero(3);
  se_ref << 0.107455, 0.107455, 0.107455;
  orb.BSESinglets().eigenvalues() = se_ref;
  // reference coefficients
  Eigen::MatrixXd spsi_ref = Eigen::MatrixXd::Zero(60, 3);
  spsi_ref << -0.00178316, -0.0558332, 0.0151767, 0.00568291, 0.0149417,
      0.0556358, 0.05758, -0.00320371, -0.00502105, 0.00379328, -0.00232255,
      -0.00817518, -0.00848959, -0.000618633, -0.00376334, -0.000395809,
      -0.00899117, 0.0023708, 7.13955e-08, 3.47498e-08, -1.08537e-08,
      0.000420413, 0.011896, -0.00320024, -0.00288487, 0.00320821, 0.0115465,
      0.0119767, 0.000355172, 0.00289343, -2.55565e-08, 1.91684e-08,
      3.01647e-09, 6.84051e-09, 2.79592e-10, -1.35767e-09, 0.00420618, 0.092451,
      -0.0239374, 0.0036276, 0.0113951, 0.0471896, 0.0465325, -0.00398807,
      -0.00484251, 0.0051614, -0.0031325, -0.0112447, -0.010968, -0.000896935,
      -0.00488789, 0.000951266, 0.0239709, -0.00630323, 0.000419006, 0.0201651,
      -0.00573095, -0.00118124, -0.0269275, 0.00700955, -0.00345217, 0.00356488,
      0.0134361, 0.013156, 9.58532e-05, 0.00315613, -2.58268e-05, -0.00124098,
      0.000352706, -1.80679e-06, -8.71959e-05, 2.47799e-05, -0.0150073,
      0.0114186, 0.0443012, 0.0180327, -0.026287, 0.0734351, -0.0643994,
      0.0257628, 0.0132084, -0.0123339, 0.0092297, -0.0148779, 0.0122493,
      -0.00246837, -0.0125735, -0.00375172, 0.00294872, 0.0112899, 0.00648758,
      -0.0055755, -0.0191436, 0.00433063, -0.00332619, -0.0128321, 0.0111166,
      -0.00969272, 0.0189659, -0.0160119, 0.00458659, 0.0107104, -0.000399315,
      0.000343129, 0.00117813, -2.80525e-05, 2.41086e-05, 8.2778e-05,
      -0.0450479, -0.00034974, -0.015063, 0.0655838, 0.0115163, -0.022358,
      0.0220978, 0.0583236, 0.0513097, -0.0119156, 0.00490159, 0.0116429,
      -0.0132479, -0.0146227, -0.00863565, -0.0118978, -0.000282044,
      -0.00400272, 0.0199347, 0.00139057, 0.00635067, 0.0131991, 0.000163177,
      0.00441453, 0.0159091, -0.00241207, -0.0110696, 0.0123057, 0.0171503,
      0.0119626, -0.00122682, -8.55716e-05, -0.00039083, -8.62007e-05,
      -6.0128e-06, -2.746e-05, -0.0304688, -0.954049, 0.259333, 0.0971056,
      0.255313, 0.950672, 0.983887, -0.0547431, -0.0857965, 0.0170489,
      -0.0104387, -0.036743, -0.0381557, -0.00278036, -0.0169143, -0.00177889,
      -0.04041, 0.0106552, -2.23782e-07, 2.40738e-07, 1.03401e-07, -0.000182535,
      -0.00516415, 0.00138942, 0.00125201, -0.00139237, -0.00501195,
      -0.00519809, -0.000154171, -0.00125602, 4.03664e-08, -6.04796e-08,
      -4.6768e-08, -2.38233e-09, 2.31605e-09, 1.35922e-09;
  orb.BSESinglets().eigenvectors() = spsi_ref;
  orb.CalcCoupledTransition_Dipoles();

  std::vector<Eigen::Vector3d> dipoles = orb.TransitionDipoles();
  std::vector<Eigen::Vector3d> dipoles_ref;
  dipoles_ref.push_back(Eigen::Vector3d(0.110512, 0.048776, -0.0515914));
  dipoles_ref.push_back(Eigen::Vector3d(-0.13408, 0.0969472, 0.0261392));
  dipoles_ref.push_back(Eigen::Vector3d(0.0586073, 0.121606, -0.0606862));

  for (Index i = 0; i < 3; i++) {
    bool check = dipoles[i].isApprox(dipoles_ref[i], 1e-5);
    BOOST_CHECK_EQUAL(check, true);
    if (!check) {
      std::cout << "ref" << i << std::endl;
      std::cout << dipoles_ref[i].transpose() << std::endl;
      std::cout << "result" << i << std::endl;
      std::cout << dipoles[i].transpose() << std::endl;
    }
  }

  Eigen::VectorXd oscs = orb.Oscillatorstrengths();
  Eigen::VectorXd oscs_ref = Eigen::VectorXd::Zero(3);
  oscs_ref << 0.001236, 0.00201008, 0.00156925;

  bool check_oscs = oscs.isApprox(oscs_ref, 1e-5);
  BOOST_CHECK_EQUAL(check_oscs, true);
  if (!check_oscs) {
    std::cout << "result" << std::endl;
    std::cout << oscs << std::endl;
    std::cout << "ref" << std::endl;
    std::cout << oscs_ref << std::endl;
  }
}

BOOST_AUTO_TEST_SUITE_END()
